长记忆时间序列定义（时域、频域和谱密度）

• 时域

  如果时间序列\({X_t}\)的自相关系数为\(\rho_\tau\)，其中\(\tau\)为滞后期数，若自相关函数的绝对值满足以下条件： \[ lim_{n\to \infty}\sum_{\tau=-n}^n\vert\rho_{\tau}\vert\to \infty \] 则称该序列为长记忆时间序列

• 频域

  如果平稳时间序列\(X_t\)的自相关函数\(\rho_{\tau}\)随滞后期数\(\tau\)的上升而缓慢变小，即满足： \[ \rho_{\tau}\sim c\tau^{2d-1},\tau\to \infty \] 其中\(c\)为常数，\(d\lt0.5\)，称该序列为长记忆时间序列

  • 当\(d\lt0\)时，\(lim_{n\to \infty}\sum_{j=-n}^n\vert \rho_j\vert\lt\infty\),称为中记忆过程
  • 当\(d\in(0,0.5)\)时，\(lim_{n\to \infty}\sum_{j=-n}^n\vert\rho_j\vert = \infty\),称为长记忆过程

• 谱密度

  若时间序列的谱密度函数\(f(\omega)\)满足：

  • \(f(\omega)\)随频率\(\omega\to 0\)而趋于无穷
  • \(f(\omega)\)在除去至多有限个\(\omega\)的值以外的所有其他\(\omega\)的值上有界

长记忆性的检验方法:计算Hurst指数的R/S方法和通过假设检验的修正的R/S检验

以下两种方法通过求解Hurst指数的值，可以进一步通过\(H = d+\frac{1}{2}\)确定\(\hat d\).

• 1、R/S分析法（又称重标极差统计法）（使用计算Hurst指数的方法检验长记忆性）

  基本思路是将时间序列\(X_t\)划分为A个长度为n的子区间，时间序列的长度为T,将子序列记作\(D_a,a=1,2,...,A\).子序列\(D_a\)的元素为\(x_{i,a},i,2,...,n\),均值为\(\bar x_a\).计算\(D_a,a=1,2,...,A\)的累计离差为：

\[ X_{k,a}=\sum_{i=1}^k(x_{i,a}-\bar x),i=1,2,...,n \]

对于每个子序列，可以得到n个累计离差，则子序列\(D_a,a=1,2,...,A\)的极差\(R_a\)可以表示为：\(R_a = max(X_{k,a})-min(X_{k,a})\)

重标极差表示为：\(\tilde Q_n = (R/S)_n = \frac{1}{A}\sum_{a=1}^A \frac{R_a}{S_a}\)

其中，\(R_a\)和\(S_a\)分别为子区间a内部数据的极差和标准差。当n取不同的值，以下关系成立：\((R/S)_n = kn^H\)

其中H为\(Hurst\)指数,k为常数。

故此通过两边求对数得到\(ln\tilde Q_n = Hlnk\),通过回归n取不同值下产生的多组数据可以拟合出\(H\)的取值，由\(Hurst\)指数来判断时间序列的长记忆性使用下列的判别标准：

• 当\(0\le H\lt0.5\),序列具有反持久性，被称作均值回复的分数布朗运动，即时间序列存在返回起始点的倾向趋势，\(H\)的值越接近零，返回起始点的可能性越大。

• 当\(H=0.5\)时，序列是标准的随机游走，变量之间不是自相关的，可以认为序列只存在短记忆性，不存在长记忆性,为短程过程

• 当\(0.5\le H \lt 1\)时，说明时间序列存在长记忆性，并且\(H\)值越大，长记忆性越强

• 当\(H=1\)时，序列可以用现在预测未来，具有长记忆性

• 2、修正的R/S分析法（使用构建统计量的方法进行长记忆性的检验）

R/S方法在检验时间序列中的长记忆方面也表现得很好，但是学者研究发现，该方法有一些弱点，其中最重要得缺陷是，无法准备区分短记忆性和长记忆性（对时间序列当中的短程记忆比较敏感）。当序列中只有短期记忆性或同时存在短记忆性和长记忆性的时候，R/S分析法往往会做出错误的判断。为了克服\(R/S\)分析法的缺陷，提出下面修正的\(R/S\)方法。

修正的R/S统计量：\(Q_n(q) = \frac{R(n)}{\sigma_n(q)}\)

\[ \hat \sigma^2_n(q) = \frac{1}{n}\sum_{j=1}^k(X_j-\bar X_n)^2+\frac{2}{n}\sum_{j=1}^n\omega_j(q)[\sum_{t=j+1}^n(X_j-\bar X_n)(X_{t-j}-\bar X_n)] \equiv \hat\gamma_0+2\sum_{j=1}^q\omega_j(q)\hat\gamma_j \]

\[ \omega_j(q) = 1-\frac{j}{q+1} \]

由于可以计算统计量\(V_n(q) = \frac{Q_n(q)}{\sqrt{n}}\):

对于短程过程：V是定义在\([0,1]\)上的布朗桥极差，\(E(V)=\sqrt{\pi/2}\approx 1.25,std(V)=\sqrt{\pi(\pi-3)/6}\approx 0.27\)

对于长记忆过程，\(V\to \infty\)

即如果时间序列只有短自相关而没有长记忆性，那么\(V\)的极限分布是稳定的，可以通过\(V\)的极限分布情况来对长记忆性进行检验。

R/S分析的原假设是模型没有长记忆性。

• ARFIMA模型特例：ARFIMA(0,d,0)模型及其性质

如果对于过程ARFIMA(0,d,0),时间序列\(X_t\)还是下面差分方程的零均值平稳解，\(d\in(-0.5,0.5)\) \[ (1-B)^dX_t = \epsilon_t,\epsilon\sim WN(0,\sigma^2) \] 则称为分数次求和白噪声。

如果\(d\in(-0.5,0.5)\),则 \[ (1-B)^dX_t = \epsilon_t \] 存在唯一的纯非确定平稳解\(X_t\),且 \[ X_t = \sum_{k=0}^{\infty}\psi_k\epsilon_{t-k} \equiv (1-B)^{-d}\epsilon_t \] 这里： \[ \psi_k = \frac{\Gamma(k+d)}{\Gamma(d)\Gamma(k+1)}=\prod_{0\le j\le k}\frac{j-1+d}{j}\sim c\times k^{-1-d} \] 并且有：方差函数、自相关函数、偏自相关函数： \[ \gamma(0)=\frac{\sigma^2\Gamma(1-2d)}{\Gamma^2(1-d)} \]

\[ \rho_k = \frac{\Gamma(k+d)\Gamma(1-d)}{\Gamma(k-d)\Gamma(d)}=\prod_{0\le j\le k}\frac{j-1+d}{j} \]

\[ \alpha(k)=\frac{d}{k-d},k=1,2,... \]

由\(\Gamma(x)\sim\sqrt{2\pi}e^{-x+1}(x-1)^{x-0.5}\)得\(\rho_k\sim \frac{k^{2d-1}\Gamma(1-d)}{\Gamma(d)}\)

当\(0\lt d\lt 0.5\)时，分数次求和白噪声是长记忆过程(\(\rho_k\to \infty\))。

• 引入滞后算子表示\(ARFIMA(p,d,q)\)过程及其性质

\(ARFIMA\)过程可以写成：\(\phi(B)(1-B)^dX_t = \theta(B)\epsilon_t,\epsilon\sim i.i.d WN(0,\sigma^2),d\in(-0.5,0.5)\),其中\((1-B)^d\)称为分数次差分算子.从模型可知，\(ARFIMA\)模型同时考虑序列中的长期和短期记忆性问题，这个模型使用\(p+q\)个参数描述序列中的短期记忆性，利用参数\(d\)描述序列中的长期记忆性。

性质：

• 如果\(X_t\)是\(ARFIMA(p,d,q)\)过程，\(d\lt0.5\),那么，\((1-B)X_t\)是\(ARFIMA(p,d-1,q)\)过程，由于若\(d\in(0,0.5)\),则\(d-1\in(-1,-0.5)\),故此，算子\((1-B)\)可以将长记忆过程变为中记忆过程。

• 当\(d\le -0.5\)时，有唯一平稳解，但是此解不可逆

• 当\(-0.5\lt d\lt0\)时，有唯一平稳解，且可逆，中记忆过程

• 当\(d=0\)时，有唯一平稳解，且可逆，短记忆过程

• 当\(0\lt d\lt 0.5\)时，有唯一平稳解，且可逆，长记忆过程

• 当\(d\ge0.5\)时，其解非平稳

• 由上可得，当\(\vert d\vert \lt0.5\)的时候，\(ARFIMA\)过程平稳可逆，进一步有以下结论：当\(-0.5\lt d\lt 0.5,d\neq0\)时，自相关函数为： \[ \rho_k = c\times k^{2d-1},d\neq0 \] 自相关函数呈现双曲律递减。

对于平稳可逆的\(ARFIMA(p,d,q)\),可以进一步写成：\(\phi(B)(1-B)^dX_t = \theta(B)\epsilon_t\),满足：

1、平稳过程的充要条件是：\(\phi(z)=0\)的根全在单位圆外

2、可逆过程的充要条件是：\(\theta(z)=0\)的根全在单位圆外

3、\(\rho_k = c\times k^{2d-1},k\to \infty,c\neq0\)

• \(ARFIMA\)过程的参数估计：均值估计，自相关系数估计和分数次\(d\)和系数参数的估计（见\(PPT\)）

• \(ARFIMA\)模型的预测：使用建立后的模型进行预测

• \(ARFIMA\)模型的建模过程

  • 1、数据预处理并进行正态性检验，清除原始数据中的趋势和波动影响

  • 2、滤除序列中的短记忆因素，为了突出序列的长记忆性，需要将短记忆因素滤除，可以通过\(AR\)模型来实现

  • 3、分数阶差分，通过估计参数\(d\)的方法（如\(R/S\)方法）实现参数估计，并完成差分得到差分后的序列

  • 4、接下来的步骤就变成一个普通的\(ARMA\)模型参数估计，利用第三步得到的新序列进行p，q定阶，并利用最小二乘估计或者是其他方法确定模型中的\(p+q\)个参数

1.判断时间序列是否长记忆

1.2 重新标度极差统计量(rescaled-range)方法

\[Q_T = R_T/S_T\]

• 短期关联过程：\(Q_T=O_p(T^{1/2})\)

• 长记忆过程： \(Q_T=O_p(T^{H})\)

其中H=d+1/2称为Hurst指数。 在 \(ln(Q_T)\) 对 \(ln(T)\) 的散点图上，短期记忆 过程的点应分布在斜率1/2的直线附近， 长记忆过程的点对应的直线斜率大于1/2. 回归方法得到对Hurst指数的估计(即斜率估计，d的估计为Hurst指数估计-0.5）。

实际中用修正的R/S分析进行检验。

计算R/S的python代码:

import numpy as np
import pandas as pd
def hurst(ts):
    ts = list(ts)
    N = len(ts)
    if N < 20:
        raise ValueError("Time series is too short! input series ought to have at least 20 samples!")

    max_k = int(np.floor(N/2))
    R_S_dict = []
    for k in range(10,max_k+1):
        R,S = 0,0
        # split ts into subsets
        subset_list = [ts[i:i+k] for i in range(0,N,k)]
        if np.mod(N,k)>0:
            subset_list.pop()
            #tail = subset_list.pop()
            #subset_list[-1].extend(tail)
        # calc mean of every subset
        mean_list=[np.mean(x) for x in subset_list]
        for i in range(len(subset_list)):
            cumsum_list = pd.Series(subset_list[i]-mean_list[i]).cumsum()
            R += max(cumsum_list)-min(cumsum_list)
            S += np.std(subset_list[i])
        R_S_dict.append({"R":R/len(subset_list),"S":S/len(subset_list),"n":k})
    
    log_R_S = []
    log_n = []
    print(R_S_dict)
    for i in range(len(R_S_dict)):
        R_S = (R_S_dict[i]["R"]+np.spacing(1)) / (R_S_dict[i]["S"]+np.spacing(1))
        log_R_S.append(np.log(R_S))
        log_n.append(np.log(R_S_dict[i]["n"]))

    Hurst_exponent = np.polyfit(log_n,log_R_S,1)[0]
    return Hurst_exponent

• 计算R/S的R代码

Hurst Exponent: Calculates the Hurst exponent using R/S analysis.pracma包中有hurstexp函数。

2.模型形式

library(readr)
library(xts)
library(forecast)
set.seed(1234)
sim <- fracdiff::fracdiff.sim(n=1500, d= 0.3)$series
ts.plot(sim)

par(mfrow=c(1,2))
acf(sim)
pacf(sim)

sim <- fracdiff::fracdiff.sim(n=1500, d=0.2,ar=0.7,ma=-0.5)$series
ggtsdisplay(sim)

ARFIMA建模示例

da <- read_table2(
  "~/Downloads/ftsdata/d-ibm3dx7008.txt",
  col_types=cols(.default=col_double(),
                 Date=col_date("%Y%m%d"))
)
xts.crspw <- xts(da[,-1], da$Date)
vw <- abs(da$vwretd)
ew <- abs(da$ewretd)
rm(da)
acf(vw, main="", lag.max=300)

acf(ew, main="", lag.max=300)

round(mean(vw),1)

## [1] 0

#极大似然估计
mres <- fracdiff::fracdiff(vw, nar=1, nma=1)
summary(mres)

## 
## Call:
##   fracdiff::fracdiff(x = vw, nar = 1, nma = 1) 
## 
## Coefficients:
##    Estimate Std. Error z value Pr(>|z|)    
## d  0.490938   0.007997   61.39   <2e-16 ***
## ar 0.113389   0.005988   18.94   <2e-16 ***
## ma 0.575895   0.005946   96.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma[eps] = 0.0065619 
## [d.tol = 0.0001221, M = 100, h = 0.0003742]
## Log likelihood: 3.551e+04 ==> AIC = -71021.02 [4 deg.freedom]

fracdiff::fdGPH(vw)

## $d
## [1] 0.372226
## 
## $sd.as
## [1] 0.0698385
## 
## $sd.reg
## [1] 0.06868857

residuals <- residuals(mres)
Box.test(residuals,lag=10,type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  residuals
## X-squared = 42.354, df = 10, p-value = 6.485e-06

Box.test(residuals,lag=20,type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  residuals
## X-squared = 78.826, df = 20, p-value = 6.207e-09

Box.test(residuals,lag=50,type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  residuals
## X-squared = 157.48, df = 50, p-value = 4.72e-13

Box.test(residuals,lag=100,type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  residuals
## X-squared = 235.43, df = 100, p-value = 6.211e-13

时间序列模拟

• 1. 短记忆(d=0)

set.seed(12)
sim <- fracdiff::fracdiff.sim(n=2000, d=0)$series
ggtsdisplay(sim)

• 2. 中记忆(-0.5 < d < 0)

当时，有唯一平稳解，且可逆，中记忆过程，谱密度函数在0点等于0

set.seed(12)
sim <- fracdiff::fracdiff.sim(n=2000, d=-0.2)$series
ggtsdisplay(sim)

注意：fracdiff.sim取值为[-0.5,0.5],否则会报错。

• 3. d>=0.5,解非平稳，d<=-0.5,解不可逆

da <- read_table2("~/Downloads/ftsdata/m-unrate.txt", col_types = cols(.default = col_double()))
xts.unrate <- xts(
  da[["rate"]], with(da, make_date(Year, mon, dd)))
indexClass(xts.unrate) <- "yearmon"
ts.unrate <- ts(da[["rate"]], start=c(1948,1), frequency = 12)
ts.dur <- diff(ts.unrate)
#作时序图
ts.unrate %>% autoplot()

ts.dur %>% autoplot()

ggAcf(xts.unrate,lag.max=300)

ggPacf(xts.unrate,lag.max = 300)

#acf和pacf中体现出其长期趋势，比较大的lag,
#考虑用arfima建模

round(mean(xts.unrate),1)

## [1] 5.7

#极大似然估计
xts.unrate <- xts.unrate[!is.na(xts.unrate)]
mres <- fracdiff::fracdiff(xts.unrate[-c(1:4)], nar=3, nma=3)
summary(mres)

## 
## Call:
##   fracdiff::fracdiff(x = xts.unrate[-c(1:4)], nar = 3, nma = 3) 
## 
## Coefficients:
##      Estimate Std. Error z value Pr(>|z|)    
## d    0.076760   0.001276  60.158   <2e-16 ***
## ar1  0.867848   0.032723  26.521   <2e-16 ***
## ar2  0.903783   0.034787  25.981   <2e-16 ***
## ar3 -0.789111   0.038956 -20.257   <2e-16 ***
## ma1 -0.061286   0.030886  -1.984   0.0472 *  
## ma2  0.637042   0.040849  15.595   <2e-16 ***
## ma3 -0.279579   0.026591 -10.514   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma[eps] = 0.1958061 
## [d.tol = 0.0001221, M = 100, h = 1.681e-06]
## Log likelihood: 158.6 ==> AIC = -301.1344 [8 deg.freedom]

residuals <- residuals(mres,standardize=T)
checkresiduals(residuals)

Box.test(residuals,lag=10,type = "Ljung")#残差白噪声

## 
##  Box-Ljung test
## 
## data:  residuals
## X-squared = 10.101, df = 10, p-value = 0.4316

Box.test(residuals,lag=15,type = "Ljung")#残差白噪声

## 
##  Box-Ljung test
## 
## data:  residuals
## X-squared = 21.716, df = 15, p-value = 0.1155

Box.test(residuals,lag=20,type = "Ljung")#残差非白噪声 模型已经ok了

## 
##  Box-Ljung test
## 
## data:  residuals
## X-squared = 32.574, df = 20, p-value = 0.03756

Box.test(residuals^2,lag=10) #残差平方非白噪声

## 
##  Box-Pierce test
## 
## data:  residuals^2
## X-squared = 102.03, df = 10, p-value < 2.2e-16

library(fGarch)
mod1 <- garchFit( ~ garch(1,1), data= residuals, trace=FALSE)
mod1

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = residuals, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x7ffcea5b8a98>
##  [data = residuals]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu        omega       alpha1        beta1  
## -0.00216283   0.00072347   0.07851152   0.90367406  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -0.0021628   0.0061894   -0.349  0.72676    
## omega   0.0007235   0.0004522    1.600  0.10963    
## alpha1  0.0785115   0.0239544    3.278  0.00105 ** 
## beta1   0.9036741   0.0313072   28.865  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  213.6704    normalized:  0.2852742 
## 
## Description:
##  Mon Mar  8 21:12:09 2021 by user:

#其中1表示均值方程是一个常数，
#输出结果中mu表示均值方程的均值，omega表示alpha0,alpha1为alpha1
resi <- residuals(mod1, standardize=TRUE)
checkresiduals(resi)

#标准化残差的白噪声检验
Box.test(resi,lag=12)

## 
##  Box-Pierce test
## 
## data:  resi
## X-squared = 17.621, df = 12, p-value = 0.1277

Box.test(resi^2,lag=12)

## 
##  Box-Pierce test
## 
## data:  resi^2
## X-squared = 18.377, df = 12, p-value = 0.1047

ggAcf(resi)

ggPacf(resi)

#经过garch模型以后，残差白噪声，模型ok.